# Replication Archive for: 
# Coppock, Alexander and Donald P. Green. 2020. 
# "Do Belief Systems Exhibit Dynamic Constraint?" 
# The Journal of Politics, Forthcoming.

library(tidyverse)
library(estimatr)
library(coefplot)

source("programs/helpers/helper_functions.R")
study_1_mturk <- read_rds("data/clean/study_1_mturk_cleaned.rds")
study_1_elite <- read_rds("data/clean/study_1_elite_cleaned.rds")
study_1_lucid <- read_rds("data/clean/studies_1_2_replications_lucid_cleaned.rds")

study_1_mturk <- study_1_mturk %>% filter(Z != "climate")
study_1_lucid <- study_1_lucid %>% filter(Z != "climate")

dv_df <- read_rds("data/raw/dv_df.rds")
dv_df_mturk <- filter(dv_df, dv_cat != "climate")
dv_df_elite <- filter(dv_df, dv_cat != "climate")
dv_df_lucid <- filter(dv_df, dv_cat != "climate")


# regressions -------------------------------------------------------------

all_fits_mturk <-
  dv_df_mturk$dvs %>%
  map( ~ formula(paste0(., "_w1 ~ Z"))) %>%
  map( ~ lm_robust(., data = study_1_mturk)) %>%
  map(tidy) %>%
  bind_rows 


all_fits_elite <-
  dv_df_elite$dvs %>%
  map( ~ formula(paste0(., "_w1 ~ Z"))) %>%
  map( ~ lm_robust(., data = study_1_elite)) %>%
  map(tidy) %>%
  bind_rows 

all_fits_lucid <-
  dv_df_lucid$dvs %>%
  map( ~ formula(paste0(., "_w1 ~ Z"))) %>%
  map( ~ lm_robust(., data = study_1_lucid)) %>%
  map(tidy) %>%
  bind_rows 



# prepare for plotting ----------------------------------------------------

all_fits <-
  bind_rows(`Elite Sample` = all_fits_elite,
            `MTurk Sample` = all_fits_mturk,
            `Lucid Sample` = all_fits_lucid,
            .id = "Study") %>%
  mutate(dvs = gsub(outcome, pattern = "_w1", replacement = ""),
         Study = factor(Study, levels = c("MTurk Sample",
                                  "Elite Sample",
                                  "Lucid Sample")))

gg_df <-
  all_fits %>%
  filter(term != "(Intercept)") %>%
  left_join(dv_df) %>%
  mutate(treatment = factor(term,
                            levels = c("Zpaul", "Zwallstreet", "Zveterans", "Zamtrak"),
                            labels = c("Treatment: Flat Tax", "Treatment: Wall Street", "Treatment: Veterans", "Treatment: Amtrak")),
         dv_cat_label = factor(dv_cat,
                               levels = c("flat", "wall", "vets","amtrak"),
                               labels = c("Outcomes: Flat Tax", "Outcomes: Wall Street", "Outcomes: Veterans", "Outcomes: Amtrak")),
         dv_labels_rev = factor(dv_labels, levels = rev(levels(dv_labels))))


highlight_df <-
  tibble(
    treatment = c(
      "Treatment: Wall Street",
      "Treatment: Amtrak",
      "Treatment: Veterans",
      "Treatment: Flat Tax"
    ),
    dv_cat_label = c(
      "Outcomes: Wall Street",
      "Outcomes: Amtrak",
      "Outcomes: Veterans",
      "Outcomes: Flat Tax"
    )
  )

g <-
  ggplot(gg_df,
         aes(
           x = estimate,
           y = dv_labels_rev,
           group = Study,
           color = Study,
           shape = Study
         )) +
  geom_rect(
    data = highlight_df,
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = -Inf,
      ymax = Inf,
      x = NULL,
      y = NULL,
      group = NULL,
      color = NULL,
      shape = NULL
    ),
    fill = 'grey',
    alpha = 0.4,
    show.legend = FALSE
  ) +
  geom_point(position = position_dodgev(height = -.5)) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 position = position_dodgev(height = -.5),
                 height = 0) +
  scale_color_manual(values = bpr_colors) +
  scale_x_continuous(breaks = c(0, 0.33, 0.67)) +
  facet_grid(dv_cat_label ~ treatment, scales = "free_y") +
  coord_cartesian(xlim = c(-.2, .8)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  theme(
    axis.title = element_blank(),
    legend.position = "bottom",
    strip.background = element_blank(),
    legend.key.width = unit(4, "lines"),
    text = element_text(size = 8)
  )
g

ggsave("../Drafts/figures/study_1.pdf", g, height = 5, width = 6.5)
ggsave("../Drafts/figures/study_1.png", g, height = 5, width = 6.5)



# in-text figures re: significance -----------------------------------------

gg_df <-
  gg_df %>%
  mutate(
    dv_cat_harmony = str_split_fixed(pattern = ": ", as.character(dv_cat_label), n = 2)[, 2],
    treatment_harmony = str_split_fixed(pattern = ": ", as.character(treatment), n = 2)[, 2],
    same_cat = if_else(dv_cat_harmony == treatment_harmony, "Within Domain", "Across Domain")
  )


gg_df %>%
  group_by(same_cat) %>%
  mutate(
    p_bon = p.adjust(p.value, method = "bonferroni"),
    p_holm = p.adjust(p.value, method = "holm"),
    p_BH = p.adjust(p.value, method = "BH")
  ) %>%
  summarize(sum(p_bon < 0.05),
            sum(p_holm < 0.05),
            sum(p_BH < 0.05),
            sum(p.value < 0.05),
            n = n())

gg_df %>%
  group_by(Study, same_cat) %>%
  summarize(weighted.mean(estimate, w = (1/std.error^2)))

